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1 . 0  SUMMARY 

This  technical  report  presents  the  results  of  the  second 
six  months  of  research  on  numerical  methods  for  two-dimensional 
process  modeling  under  Contract  MDA903-80-C-0498 ,  DARPA 
Order  No.  3984.  During  the  first  six  months  of  the  contract 
effort,  a  fast  algorithm  for  computing  dopant  spread  was 
found  and  tested.  In  the  current  period,  a  code  incorporating 
nonuniform  motion  of  the  oxide-silicon  boundary  was  written 
around  this  algorithm,  and  tests  were  made  on  typical  process 
conditions  which  showed  that  the  most  difficult  process  steps, 
involving  simultaneous  nonuniform  oxide  growth  and  nonlinear 
dopant  diffusion,  could  be  accurately  simulated  in  under 
30  seconds  of  CPU  time  on  the  Cyber  176.  Transfer  of  the 
basic  algorithm  to  the  process  and  device  modeling  group  at 
Stanford  has  been  initiated.  In  addition,  the  capabilities 
of  this  algorithm  were  reviewed  at  the  Second  International 
Conference  on  Numerical  Analysis  for  Semiconductor  Devices 
and  Integrated  Circuits,  held  in  Dublin,  Ireland, 

June  17-19,  1981. 

1 . 1  TASK  OBJECTIVES 

The  overall  objective  of  this  program  is  to  develop 
fast  and  accurate  methods  for  computer  modeling  of  the 
two-dimensional  spread  of  dopants  and  other  defects  during 
VLSI  circuit  fabrication.  Our  initial  goals  have  been  to 
demonstrate  these  methods  for  nonlinear  diffusion  of  a 
single  dopant  during  oxide  growth,  and  to  provide  the 
resulting  computational  algorithms  in  a  form  suitable  for 


1 


Rockwell  International 

SC5271.6SA 


incorporation  into  a  general  process  modeling  computer  code, 
such  as  Stanford's  program  SUPREM. 

Three  tasks  were  defined  for  the  first  year's  effort: 

Task  1  —  Analysis  of  Numerical  Methods 

Task  2  -  Formulation  of  Test  Cases 

Task  3  -  Algorithm  Development  and  Evaluation  . 

Their  objectives  are,  respectively: 

Task  1  —  Select  promising  algorithms  from  other  disciplines, 
primarily  fluid  dynamics,  where  considerable 
progress  on  numerical  methods  for  multidimen¬ 
sional  problems  has  recently  been  made. 

Task  2  -  a)  Determine  realistic  parameter  ranges  to 

be  used  as  test  cases  for  source  and  drain 
diffusion  into  a  short  MOSFET  channel  region. 

b)  Develop  and  analyze  two-dimensional  diffusion 
problems  for  which  dopant  profiles  can  be 
obtained  to  high  accuracy  by  existing 
techniques . 

Task  3  —  a)  Evaluate  each  selected  algorithm  for  speed 
and  accuracy. 

b)  Adapt  the  most  promising  algorithm (s)  to 
solve  the  combined  oxidation-nonlinear 
diffusion  problem  for  a  single  dopant 
species  in  two  dimensions. 

All  objectives  of  the  first  year's  effort  have  been  accomplished. 
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The  specific  objectives  for  the  second  year  include: 

1.  Effective  transfer  of  the  basic  algorithm  to 
the  integrated  circuits  community; 

2.  Extension  of  the  code  to  treat  multiple 
interacting  species  and  three-dimensional 
redistribution;  and 

3.  Exploration  of  the  computational  requirements 
posed  by  better  physical  models  for  the 
underlying  processes  of  chemical  reaction 
and  defect  generation  and  migration. 

Much  of  the  groundwork  for  accomplishing  the  first  two 
objectives  has  already  been  laid.  A  computer  code  incorporating 
the  basic  algorithm  and  providing  user  interface  routines 
adapted  from  SUP REM  is  in  the  final  stages  of  documentation, 
and  is  presently  scheduled  for  delivery  to  Stanford  in  August. 

It  will  also  be  available  directly  to  interested  users. 

With  regard  to  multiple  species  and  three-dimensional 
redistribution,  only  the  addition  of  subroutines  to  the 
existing  code  will  be  required,  as  the  basic  algorithm  is 
already  set  up  to  handle  the  more  general  cases. 

1.2  TECHNICAL  PROBLEM 

The  fabrication  of  VLSI  devices  requires  production  of 
features  of  submicron  size  and  separation.  Electrical 
characteristics  such  as  threshold  and  punchthrough  voltages 
will  be  sensitive  to  dopant  spread  into  critical  areas 
adjacent  to  the  original  features.  Experimental  control 
of  this  spreading,  without  guidance  from  accurate  computer 
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modeling,  will  be  costly,  tedious,  and  time-consuming. 
However,  the  use  of  standard  numerical  methods  to  achieve 
an  adequate  modeling  capability  is  also  costly  and  time- 
consuming.  One  should  therefore  seek  advanced  methods, 
drawn  from  areas  such  as  fluid  dynamics,  where  considerable 
effort  and  ingenuity  have  been  expended  in  recent  years 
to  develop  fast  and  accurate  solvers  for  the  characterization 
of  multidimensional,  time-dependent  phenomena. 

1.3  GENERAL  METHODOLOGY 

Based  upon  our  own  ongoing  research  in  computational 
nonlinear  aerodynamics,  we  identified  several  promising 
approaches  to  the  development  of  a  fast  solver  for  two- 
dimensional  diffusion  problems.  After  a  preliminary 
screening,  a  few  of  these  were  selected  for  adaptation  to 
the  problem  of  dopant  spread  during  oxidation  or  annealing. 
These  algorithms  were  tested  for  speed  and  accuracy  on  the 
problem  of  nonlinear  dopant  diffusion  into  the  channel 
region  of  a  MOSFET,  as  well  as  on  simpler  problems  for 
which  the  actual  dopant  profiles  could  be  accurately 
obtained  by  other  means . 

1.4  TECHNICAL  RESULTS 

The  principal  technical  result  obtained  during  the 
second  six  months  of  contract  effort  was  the  successful 
incorporation  of  nonuniform  oxide  growth  into  the  dopant 
redistribution  code.  Typical  process  steps  involving 
implant  redistribution  in  the  vicinity  of  a  growing  bird's 
beak  field  oxide  have  been  run  in  under  30  seconds  on  the 
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Cyber  176  (see  Appendix) .  Given  computer  times  of  this 
order,  the  iterative  use  of  this  code  to  design  a  fabrication 
process  in  two  dimensions,  much  as  SUPREM  is  used  in  one 
dimension,  now  appears  practical. 

Additional  results  worthy  of  note  are  the  verification 
of  accuracy  for  the  method  of  lines  algorithm,  reported  in 
Section  2  below,  and  a  comparison  of  those  factors  which 
dominate  the  time  and  storage  requirements  for  this  algorithm, 
and  for  the  finite-element  algorithm  now  in  common  use  for 
solving  the  partial  differential  equations  encountered  in 
process  and  device  simulation,  given  in  the  Appendix. 

1.5  IMPORTANT  FINDINGS  AND  CONCLUSIONS 

The  first  year's  effort  has  produced  a  fast  and  accurate 
algorithm  which  handles  an  essential  part  of  the  VLSI  process 
simulation  problem:  dopant  spread  in  two  dimensions  during 
annealing  or  nonuniformly  oxidizing  process  steps.  Its 
speed  is  sufficiently  high  to  encourage  iterative  use  by  the 
process  engineer  in  selecting  process  parameters  for  those 
steps.  While  the  code  should  be  useful  in  its  present  form, 
which  includes  user  interface  routines  adapted  from  SUPREM, 
it  is  intended  to  function  within  the  framework  of  a  complete, 
two-dimensional  process  simulator.  By  providing  this  code 
to  the  Stanford  modeling  group,  we  hope  to  stimulate  the 
creation  of  such  a  simulator.  Alternatively,  process  design 
groups  presently  engaged  in  VLSI  fabrication  can  utilize 
the  code  as  it  stands. 

1.6  SPECIAL  COMMENTS 

As  part  of  the  strategy  to  achieve  effective  transfer 
of  our  results  to  the  integrated  circuits  community,  we  are 
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publicizing  the  capabilities  of  the  algorithm.  A  short 
paper  on  this  subject  was  presented  in  June  at  the  Second 
International  Conference  on  Numerical  Analysis  of  Semicon¬ 
ductor  Devices  and  Integrated  Circuits,  and  another  is 
scheduled  for  the  1981  Symposium  on  VLSI  Technology  this 
September.  Copies  of  these  papers  are  included  in  the 
Appendix  to  this  report.  Judging  from  the  response  to 
our  first  presentation,  there  will  be  a  number  of  groups 
requesting  copies  of  the  code  and  the  associated  documenta¬ 
tion  in  the  near  future.  Representatives  from  Honeywell, 
IBM,  Phillips,  Thomson-EFCIS ,  and  others  have  already 
expressed  interest. 

1.7  IMPLICATIONS  FOR  FURTHER  RESEARCH 

One  essential  advantage  of  the  basic  algorithm  can  be 
stated  as  follows:  Both  time  and  storage  requirements  grow 
only  as  the  number  of  unknowns  in  the  problem.  This  implies 
that  introducing  a  second  dopant  will,  for  a  fixed  mesh, 
only  double  the  CPU  time,  or  that  treating  dopant  spread 
into  a  third  dimension  will  multiply  this  time  by  the  number 
of  intervals  chosen  along  this  dimension.  It  thus  appears 
that  the  planned  extensions  to  these  more  complex  process 
conditions  will  not  place  excessive  demands  on  the  computer. 

The  growing  availability  of  parallel-architectured 
computers  such  as  the  CDC  205  and  the  Cray  series  presents 
the  interesting  possibility  of  achieving  a  large  increase 
in  the  speed  of  the  basic  algorithm.  The  iterative  solution 
process  which  dominates  the  computer  time  can  be  recast  to 
take  advantage  of  both  pipeline  structure  and  simultaneous, 
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independent  arithmetic  operations.  Such  a  gain  in  speed 
could  bring  the  investigation  of  detailed  multiple  defect 
models  within  the  scope  of  routine  computer  analysis  on 
parallel  machines. 
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2.0  CORRELATION  WITH  ANALYTICAL  PREDICTIONS 


In  this  section,  numerical  simulations  for  the  thermal 
redistribution  of  40  keV  boron  and  arsenic  implants  are 
correlated  with  analytical  predictions.  These  simulations 
were  carried  out  for  low  dose  ( 5  *  10 11  cm-2)  implants  so  the 
governing  diffusion  equation  is  linear.  Also,  the  problem 
was  further  simplified  by  utilizing  a  non-oxidizing  ambient 
to  drive  the  implants  into  the  bulk  of  the  silicon  medium. 
For  drive-in  problems,  the  simplification  of  the  boundary 
value  problem  arises  from  the  fact  that  the  boundaries  are 
stationary.  Figure  1  is  a  qualitative  two-dimensional 
sketch  for  defining  the  drive-in  boundary  value  problem. 

In  this  figure  a  window,  y  <  |a|,  has  been  etched  out  of 
an  impenetrable  mask  along  the  z  direction.  Through  this 
window,  an  implant  is  introduced  into  the  silicon  medium 
via  a  thin  (v  0.005  ym)  protective  oxide  (SiC^)  pad.  On 
entering  the  silicon  medium,  the  implant  will  be  distributed 
two  dimensionally  as  follows:^- 


N .  x  io“ 

N(x,y,0)  =  - — -  exp 

2  \J  2 tt  Cp 


V2 

2°P 


1  1 

( 

- 

.  I 

r 

(y  +  a) 

/2  oL_ 

-  erf 


(y  -  a) 

/?  oL 


( cm" 3 )  , 


(1) 


where  N,  is  the  dose  in  cm-2,  R  -  R  -U  ,  U  is  the  oxide 
d  p  p  o  o 

pad  thickness,  R  is  the  projected  range,  and  a  and  oT  are 
P  P  ■*-' 

the  vertical  and  lateral  standard  deviations,  respectively. 
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When  the  window  is  brought  into  contact  with  a  non-oxidizing 
ambient,  the  implant  distribution  given  by  Eq.  (1)  will  be 
thermally  driven-in.  The  governing  equation  for  the  redis¬ 
tribution  of  impurities  is  obtained  on  combining  the  continuity 

2 

equation  with  Fick's  law  for  the  impurity  flux  density 
vector.  This  gives  rise  to  the  result, 

h  [D  h N]  +  &  [D  h N]  '  it N  - 0  -  121 

where  D  is  the  diffusion  coefficient,  which  is  constant  for 
low  dose  implants  and  concentration  dependent  for  high  dose 
implants.  On  the  boundaries  of  the  domain  |o  <  x  <  Lq,  y  <  bj 
of  Fig.  1  on  which  Eq.  (2)  is  defined,  the  impurity  concen¬ 
tration  satisfies  the  following  conditions: 

N  =  °  ,  (3) 

3X  x=0 


(4) 


and 


(5) 


The  silicon-silicon  dioxide  and  silicon-sapphire  interfaces 
located  at  x=0  and  x=Lq,  respectively,  are  assumed  to  be 
impenetrable  barriers  to  out  diffusion,  hence,  Eqs.  (3)  and 
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(4)  are  strictly  reflecting  boundary  conditions.  The 
boundary  conditions  at  y=±b  and  given  by  Eq.  (5)  imply  that 
the  structure  shown  in  Fig.  1  repeats  itself  periodically 
in  the  lateral  direction. 

The  solution  to  the  boundary  value  problem  defined 
by  Eqs.  (2)  to  (5)  for  a  specified  initial  distribution 
as  given  by  Eq.  (1)  can  be  obtained  analytically  for  a 
constant  diffusion  coefficient,  D,  corresponding  to  low 
dose  implants.  This  solution  will  be  in  the  form  of  a 
double  infinite  series  of  damped  two-dimensional  trigonometric 
modes  whose  generalized  Fourier  coefficients  must  be  evaluated 
numerically  in  terms  of  the  specified  initial  distribution. 
Since  this  solution  requires  a  considerable  amount  of 
programming  effort  in  order  to  simulate  the  redistribution 
of  drive-in  impurity  profiles  on  a  computer,  a  simplification 
of  the  boundary  value  problem  was  considered.  The  simplified 
problem  retains  the  boundary  condition  given  by  Eq.  (3) ,  but 
considers  the  silicon-sapphire  interface  and  the  lateral 
positions  at  y=±b  to  be  remotely  located  so  that  no  impurities 
are  reflected  from  these  boundaries.  An  excellent  approxi¬ 
mation  to  this  latter  problem  involves  the  solution  to 
Eq.  (2)  as  an  initial  value  problem  with  N(x,y,0)  specified 
by  Eq.  (1),  such  that  N(x,y,t)  satisfies  Eq.  (3)  for  all 
values  of  time  and  natural  boundary  conditions,  that  is, 
both  the  impurity  concentration  and  its  derivatives  vanish 

as  the  boundaries  located  at  x=L  and  y=±b  recede  to 

o 

infinity.  The  solution  to  this  problem  is  formally  given  as 

/•+<»  /•+<*> 

N (x,y , t)  =  J  dyoJ  dxQN (xQ,yo, 0) G (xQ,yo, 0 | x,y , t)  (6) 

*»co  0 
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where  the  kernel  under  the  integral,  given  explicitly  by 


G(xo'yo'0lx,y,t) 
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is  the  Green's  function  or  impulse  response;  that  is,  it  is 
the  response  at  the  field  point,  (x,y) ,  for  a  time  t  >  0  due 
to  an  impulse  of  impurities  which  is  introduced  at  the  source 
point,  (xo,yQ)  ,  at  time  t=0. 

The  introduction  of  Eqs.  (1)  and  (7)  into  Eq.  (6)  for 
N(xo,yQ,0)  and  G (x, y , t | xq , yQ , 0) ,  respectively,  yields 


N (x,y, t) 
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L(x)  and  H(y)  were  evaluated  numerically  using  the 
Gauss-Laguerre  and  Gauss-Hermite  quadratures,  respectively. 

2 . 1  BORON  IMPLANT 

The  initial  spatial  distribution  of  the  boron  implant 
as  specified  by  Eq.  (1),  but  with  a  uniform  background  of 
Nfa  =  1015cm-3  added  to  it,  is  plotted  in  Fig.  2.  In  carry¬ 
ing  out  the  calculation  based  on  Eq.  (1) ,  a  window  opening 
of  2a  =  1  vim  was  assumed  and  the  extreme  lateral  dimension 
of  the  domain,  b,  was  set  equal  to  1.5  vim.  The  silicon  film 

thickness,  L  ,  was  chosen  to  be  1  pm.  For  a  40  keV  boron 
o 

implant,  the  projected  range  statistics  of  Ref.  3  yielded 

the  following  parameter  values:  R  =  0.1302  pm,  a  =  0.0443  pm, 

P  P 

and  oL  =  0.0682  pm.  The  dose  of  the  boron  implant  was  low 
(5  x  10n  cm'2)  so  the  redistribution  problem  is  linear.  In 
Fig.  2(a),  the  two-dimensional  distribution  is  plotted  only 
over  the  cross-hatched  domain  of  Fig.  1  because  y=0  is  a 
plane  of  symmetry;  while  in  Fig.  2(b),  the  equi-de:»sity 
contours  are  shown.  Besides  providing  qualitative  insight 
into  the  two-dimensional  spatial  distribution  of  the  boron 
implant.  Fig.  2  likewise  provides  information  on  the  lateral 
penetration  of  the  implant. 

When  the  structure  of  Fig.  1  is  brought  in  contact  with 
an  inert  or  nonoxidizing  ambient,  such  as  argon,  the  boron 
implant  in  Fig.  2  will  be  thermally  redistributed.  For  an 
ambient  temperature  of  1000°C  the  diffusion  coefficient,  D, 
was  set  equal  to  the  non-enhanced  intrinsic  value, 

0.0055  pm2/hr,  for  boron  in  silicon.  The  redistributed 
two-dimensional  profile  corresponding  to  a  drive-in  time 
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Fig.  2(a).  Two-dimensional  Furukawa 
distribution  for  a  40  keV 
boron  implant  of  dose 

5  x  10 11  cm"2 . 
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EQUI-DENSITY  CONTOURS  (TIME  -  0.00  HRS) 
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Fig.  2(b).  The  equi-density  contours  for 
Fig.  2 (a) . 
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of  1  hr  is  shown  in  Fig.  3.  For  the  two-dimensional 
distribution  shown  in  Fig.  3(a)  and  equi-density  contours 
of  Fig.  3(b),  the  cross-hatched  computational  domain  in 
Fig.  1  was  uniformly  subdivided  into  a  grid  of  41  *  61 
points  and  the  solution  at  these  grid  points  evaluated 
by  Eq.  (8) .  One-dimensional  profiles  along  the  x  direction 
corresponding  to  y=0  (plane  of  symmetry)  and  y=a  (edge  of 
window)  are  shown  in  Fig.  4.  These  profiles  are  for  the 
situation  where  the  implant  is  decoupled  from  the  background. 
The  solid  curves  in  Fig.  4  are  the  results  predicted 
analytically  by  Eq.  (8) ,  while  the  open  and  solid  circles 
are  the  results  obtained  numerically  for  the  fine 
(Ax  =  Ay  =  0.0125  ym)  and  coarse  (Ax  =  Ay  =  0.05  ym) 
partitions,  respectively. 

An  inventory  of  the  results  shown  in  Fig.  4  reveals 
that  from  a  qualitative  viewpoint,  the  numerical  results 
for  both  the  coarse  and  fine  partitions  are  in  excellent 
agreement  with  the  analytical  results;  however,  the  numerical 
results  for  the  fine  partition  yielded  the  better  correlation. 
Although  the  numerical  results  for  the  intermediate 
(Ax  =  Ay  =  0.025  ym)  partition  are  not  plotted  in  Fig.  4, 
their  quantitative  correlation  with  the  analytical  results 
is  shown  in  Table  1.  An  inventory  of  the  results  in  Table  1 
reveals  that  analytical  and  numerical  solutions  are  in 
excellent  agreement  over  the  entire  depth  of  the  pertinent 
portion  (>  10 15  cm” 3 )  of  the  redistributed  boron  profile. 

2 . 2  ARSENIC  IMPLANT 

The  parameter  values  =  0.0269  ym,  Op  =  0.0099  ym, 
and  crT  =  0.0103  ym,  as  obtained  from  Ref.  3,  were  utilized 
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2D  DISTRIBUTION  (time=l.QO  hn) 
*=LOCJN+N.| 

R  -  IMPLANT(boron)  DISTRIBUTION 
Nfc  -  UNi  FORu'(boron)  BACKGROUND 


Y 


4  ,  0 


Fig.  3(a).  Redistributed  two-dimensional  profile  for 
the  boron  implant  shown  in  Fig.  2  after 
1  hr  of  drive-in  time  in  an  inert  argon 
ambient  at  1000°C. 
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Fig.  4.  Correlation  of  analytical  and  numerical  predictions 
along  the  x  direction  for  y=0  vm  and  y»0.5  ym  (mask 
edge)  for  the  distribution  shown  in  Fig.  3  but  with 
the  redistributed  implant  decoupled  from  the  background 


Table  1.  Correlation  of  Numerical  and  Analytical  Predictions  for  Boron  Implant* 
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♦Data:  40  keV  Boron  Implant;  Dose  5*  10ncm”z;  Drive-in  Time  =  1.0  hr;  Inert 

Ambient  (Argon)  at  1000°C. 
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with  the  same  window  opening,  that  is,  2a  =  1  ym,  to  calculate 
the  initial  distribution  for  a  low  dose  (5  x  10ncm'2)  40  keV 
arsenic  implant  according  to  Eq .  (1).  The  resultant  two- 

dimensional  distribution,  but  with  a  uniform  background 
Nb  =  1015cm-3  subtracted  from  it,  is  shown  in  Fig.  5(a), 
while  in  Fig.  5(b)  the  equi-density  contours  are  shown.  For 
a  drive-in  time  of  1  hr  in  an  inert  ambient  (argon)  at 
1000°C ,  the  redistributed  two-dimensional  arsenic  profile 
and  corresponding  equi-density  contours  are  shown  in 
Figs.  6(a)  and  6(b),  respectively.  The  results  shown  in 
Fig.  6  were  obtained  by  the  analytical  solution  of  Eq.  (8) 
using  a  non-enhanced  intrinsic  diffusion  coefficient  of 
value  0.000615  ym2/hr  for  arsenic  in  silicon.  Also,  the 
computational  domain  jo  <  x  <  Lq,  0  <  y  <  b|  was  reduced 
from  that  of  the  previous  example  by  setting  Lq  and  b  equal 
to  0.25  um  and  0.75  ym,  respectively.  The  reduction  in 
these  geometrical  lengths  was  done  such  as  not  to  impair 
the  validity  of  the  analytical  solution  given  by  Eq.  (8) . 

In  the  reduced  computational  domain,  a  fine  (Ax  =  Ay  =  0.00625  ym) 
partition  was  used  and  at  each  of  the  21 x  61  grid  points. 

Eq.  (8)  was  used  to  obtain  the  analytical  solution.  The 
results  for  the  numerical  solutions  shown  in  Fig.  7,  which 
give  the  impurity  concentration  as  a  function  of  x  for 
y  =  0  ym  and  y  =  0.5  ym  (mask  edge),  were  obtained  using 
two  rather  fine  partitions,  that  is  Ax  =  Ay  =  0.0125  ym 
and  Ax  =  Ay  =  0.00625  ym.  An  attempt  was  made  to  use  the 
intermediate  (Ax  =  Ay  =  0.025  ym)  and  coarse  (Ax  =  Ay  =  0.05  ym) 
partitions  which  were  used  previously  for  the  boron  implant. 
However,  they  yielded  numerical  results  that  departed 
considerably  from  the  analytical  results,  the  reason  being 
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♦*L0CJN-NJ 

N  -  I MPLANT(arsenic)  DISTRIBUTION 
N„  -  UNlFORM(boron)  BACKGROUND 


Fig.  5(a).  Two-dimensional  Furukawa  distribution 
for  a  40  keV  arsenic  implant  of  dose 
5  *  10 11  cm"2  . 
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EQUI-DENSITY  CONTOURS  (TIME  -  0.00  HRS) 
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Fig.  5(b).  The  equi-density  contours  for 
Fig.  5 (a) . 
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Fig.  6(a).  Redistributed  two-dimensional  profile 
for  the  arsenic  implant  shown  in 
Fig.  5  after  1  hr  of  drive-in  time  in 
an  inert  argon  ambient  at  1000°C. 
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EQUI-DENSITY  CONTOURS  (TIME  -  1.00  HRS) 
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Fig.  6(b).  The  equi-density  contours  for 
Fig.  6 (a) . 
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Fig.  7.  Correlation  of  analytical  and  numerical  predictions 
along  the  x  direction  for  y  ■  0  pm  and  y  ■  0.5  pm 
(mask  edge)  for  the  distribution  shown  in  Fig.  6, 
but  with  the  redistributed  implant  decoupled  from 
the  background . 
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that  with  these  cruder  partitions,  the  initial  two-dimensional 
distribution  of  the  low  energy  40  keV  arsenic  implant  could 
not  be  characterized  accurately  because  of  the  small  and 
cL  values. 

An  inventory  of  the  results  shown  in  Fig.  7  revealed 
that  the  numerical  results  for  both  partitions  are  in  excellent 
agreement  with  the  analytically  predicted  profiles.  The 
quantitative  correlation  between  the  analytical  and  numerical 
results,  where  the  latter  utilizes  the  Ax  =  Ay  =  0.0125  ym 
partition,  is  shown  in  Table  2. 
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Table  2.  Correlation  of  Numerical  and  Analytical  Predictions  for  Arsenic  Implant* 
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EFFICIENT  NUMERICAL  SOLUTION  OF  TWO-DIMENSIONAL  NONLINEAR 
DIFFUSION  EQUATIONS  WITH  NONUNIFORMLY  MOVING  BOUNDARIES: 

A  VERSATILE  TOOL  FOR  VLSI  PROCESS  MODELING* 

W.D.  Murphy  and  W.F.  Hall 

Rockwell  International  Science  Center,  Thousand  Oaks,  CA  91360 

C.D.  Maldonado 

Rockwell  International  Microelectronics  Research  4  Development  Center 

Anaheim,  CA  92803 

Abstract 

A  fast  algorithm*  for  the  calculation  of  two- dimensional  dopant 
redistribution  during  oxide  growth  has  been  developed.  CPU  times  on 
the  Cyber  176  are  less  than  SO  seconds  for  the  most  computationally 
difficult  process  steps ,  including  growth  of  a  bird's  beak  field 
oxide.  This  represents  a  reduction  in  computer  time  of  at  least  two 
orders  of  magnitude  over  techniques  currently  being  employed ,  such  as 
the  B-spline  finite  element  method  used  by  Varner  and  Vilson\_l'\.  Of 
greater  significance  is  the  fact  that  these  CPU  times  are  short 
enough  to  enable  practical,  iterative  computer-based  process  design 
to  be  carried  out  in  two  dimensions. 

The  lateral  diffusion  of  dopants  between  neighboring  features  on 
a  silicon  chip  can  no  longer  be  ignored  in  the  design  of  integrated 
circuit  fabrication  processes.  Because  the  size  and  separation  of 
features  in  VLSI  circuits  are  comparable  to  the  depths  over  which 
implants  are  to  be  redistributed,  their  resulting  electrical  charac¬ 
teristics  are  substantially  affected  by  the  spread  of  dopants  into 
sensitive  areas. 

In  the  Stanford  code  SUPREM,  the  LSI  process  designer  presently 
has  at  hand  a  powerful  tool  for  predicting  the  effect  of  each  process 
step  on  the  vertical  redistribution  of  dopants.  Final  dopant  pro¬ 
files  can  be  generated  in  seconds  on  a  fast  computer  such  as  the 
IBM  3033  or  Cyber  176,  enabling  the  designer  to  tailor  the  process 
quickly  and  easily  for  the  desired  result  before  expensive  and  time- 
consuming  fabrication  experiments  are  undertaken. 

• 

Our  goal  has  been  to  provide  comparable  speed  and  flexibility 
for  VLSI  process  design.  Specifically,  we  have  worked  to  develop 
fast  algorithms  for  the  calculation  of  two-dimensional  dopant  pro¬ 
files  resulting  from  VLSI  fabrication  processes,  including  nonuniform 
motion  of  the  oxide-silicon  boundary. 

As  a  measure  of  our  success,  computer  times  of  under  10  seconds 
have  been  achieved  on  the  Cyber  176  for  the  generation  of  two- 
dimensional  dopant  profiles  on  a  21*41  grid  for  a  15-minute  drive-in 


♦This  work  is  supported  in  part  by  the  Defense  Advanced  Research 
Projects  Agency  under  Contract  MDA903-80-C-0498. 
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cycle  on  boron- Implanted  silicon  at  1100°C,  with  a  peak  initial  boron 
concentration  of  1020  cm-3.  This  result  can  be  directly  compared 
with  the  work  of  Warner  and  Wilsonp],  where  similar  nonlinear  dif¬ 
fusion  problems  with  the  same  number  of  unknowns,  solved  by  second- 
order  finite  elements,  took  over  7  minutes  on  the  Cray-1  computer. 
Since  the  Cray-1  is  at  least  five  times  faster  than  the  Cyber  176, 
this  represents  an  improvement  in  speed  by  at  least  two  orders  of 
magnitude. 

Dopant  redistribution  during  high  temperature  process  cycles 
is  commonly  described  by  diffusion  of  a  single  species  with  a 
concentration-dependent  diffusivity.  Ignoring  diffusion  in  the 
oxide,  which  is  justified  in  many  cases,  one  can  characterize  the 
dopant  redistribution  problem  as  nonlinear  diffusion  in  a  bounded 
two-dimensional  domain,  with  homogeneous  boundary  conditions  at  a 
(possibly)  moving  interface: 

||=V-(DVN),  (la) 

n*DVN  =0  at  the  silicon  and  sapphire  boundaries  ,  (lb) 

n-DVN  =  n-x(k-m)(3U/3t)N  at  the  oxide  boundary  ,  (lc) 

where  U  is  the  local  oxide  thickness  (along  x),  k  is  the  segregation 
coefficient,  m  is  the  ratio  of  silicon  thickness  consumed  to  oxide 
thickness  produced,  n  is  the  unit  normal  at  each  boundary,  N  is  the 
dopant  concentration,  and  D  is  the  diffusivity. 

The  geometry  for  our  calculations  models  the  growth  of  a  bird's 
beak  field  oxide.  The  initial  condition  represents  an  ion  implant, 
which  is  modeled  by  a  Gaussian  distribution  in  depth  (x)  and  a  com¬ 
plementary  error  function  in  the  lateral  dimension  (y)  about  the  mask 
edges,  as  suggested  by  Furukawa,  et  aZ.[2].  The  rate  of  oxide  growth 
is  taken  to  vary  locally  in  such  a  way  that  there  is  a  smooth  transi¬ 
tion  from  the  high  growth  rate  on  top  of  the  field  oxide  to  the  low 
(or  zero)  growth  rate  under  the  mask.  This  is  the  general  case  for 
which  our  algorithms  are  developed;  simple  drive-in  is  obtained  by 
immobilizing  the  oxide-silicon  boundary. 

The  basic  approach  used  in  solving  the  nonlinear  diffusion  equa¬ 
tion  is  the  method  of  lines,  coupled  with  the  best  currently  avail¬ 
able  stiff  integrator,  GEARBI,  of  Hindmarsh[3].  This  integrator  is 
a  variable  order  and  variable  step  size  method  which  is  optimally 
first  to  fifth  order  accurate  in  time.  To  simplify  the  logic,  the 
changing  silicon  region  is  mapped  Into  a  rectangle  of  constant 
dimensions,  enabling  a  fixed  computational  grid  to  be  used  for  the 
discretization  and  solution  of  the  nonlinear  diffusion  equation  as 
a  system  of  ordinary  differential  equations. 

Two  types  of  mapping  have  been  used;  a  scaling  of  the  depth  which 
varies  with  lateral  position,  and  a  conformal  map  whose  level  lines 
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approximate  the  oxide-silicon  boundary  shape.  In  the  former  case, 
cross  derivatives  appear  in  the  differential  equation  and  lateral 
derivatives  appear  in  the  boundary  condition  at  the  moving  interface. 
Our  preliminary  investigations  indicate  that  relatively  fine  meshes 
and  tight  error  tolerance  are  required  to  produce  acceptable  solu¬ 
tions  for  this  approach.  In  the  latter  case,  these  tolerances  can 
be  loose,  but  time  is  consumed  in  computing  the  mapping  function  and 
its  derivatives.  In  either  case,  central  differencing  is  used  on 
the  computational  grid.  The  boundary  conditions  are  incorporated  by 
introducing  "ghost"  points  immediately  outside  the  computational 
domain  whose  values  are  determined  from  the  discretized  form  of  the 
boundary  conditions. 

The  resulting  set  of  coupled  nonlinear  ordinary  differential  equa¬ 
tions  dN/dt  =  f(N,t)  is  stiff[4],  and  consequently  an  iterative 
procedure  utilizing  the  Jacobian  3f/9N  is  required  to  converge  the 
corrector  equation  in  the  linear  multistep  method  used  to  integrate 
it.  Most  of  the  computer  time  is  used  in  solving  a  large  linear 
system  Px =  b  for  the  corrector,  where  the  matrix  P  has  the  form 
I  -  y3f/3N.  Consequently,  it  is  critical  that  efficient  methods  be 
used  in  its  solution.  In  GEARBI,  a  successive  overrelaxation 
(SOR ) [5]  technique  is  employed,  which  is  well  suited  to  our  problem 
because  3f/3N  has  only  five  nonzero  nearest-neighbor  diagonals. 

Other  advantages  of  SOR  are  minimal  storage  requirements  and,  as  P 
changes  incrementally,  fast  convergence  utilizing  the  previous 
solution  for  x  as  the  initial  iterate  for  the  next  time  step. 

CPU  time  and  storage  requirements  are  roughly  linear  with  the 
number  of  equations  being  solved.  This  is  a  great  improvement  over 
banded  matrix  techniques,  which  require  time  and  storage  proportional 
to  L(M) 3  and  to  L*1(3M+1),  respectively,  where  L  is  the  larger  and 
M  the  smaller  dimension  on  the  computational  grid.  Furthermore, 
extensions  to  coupled  species  diffusion  and  to  three  spatial 
variables  are  easily  accomplished  with  this  software. 

CPU  times  for  all  the  cases  we  have  run,  including  two-hour  boron 
drive-in  at  1000°C,  are  less  than  30  seconds  on  the  Cyber  176,  using 
computational  grids  of  up  to  31x51  points.  Results  for  the  most 
general  case  of  simultaneous  diffusion  and  nonuniform  oxide  growth 
are  illustrated  in  Figure  1.  The  conditions  for  this  case  model  an 
80  keV  boron  implant  of  3*10n  cm-2  dose,  oxidized  in  steam  at 
1000°C  for  half  an  hour.  This  problem  was  solved  on  a  31  *  51  grid 
in  12.5  seconds. 

A  typical  drive-in  problem  without  oxidation  is  illustrated  in 
Figure  2.  Here  a  high-dose  40  keV  boron  implant  is  redistributed 
over  2  hours  at  1000°C.  This  problem  required  22.5  seconds  CPU 
time  for  a  21  *  31  grid. 

Our  primary  objective  in  this  effort  is  to  provide  the  integrated 
circuits  industry  with  a  practical  means  of  extending  process  design 
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to  two  dimensions.  The  algorithm  described  above  meets  this  need. 
There  remains,  however,  the  task  of  incorporating  these  numerical 
methods  into  a  complete  process  modeling  code,  after  the  fashion  of 
SUPREM.  We  are  presently  engaged  with  the  process  modeling  group 
at  Stanford  to  achieve  this  end. 
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ABSTRACT 

The  nonlinear  and  nonuniform  moving  boundary  value  problem  which  governs 
the  redistribution  of  a  field  implant  during  the  growth  of  a  bird's  beak,  is  formu¬ 
lated  and  solved  by  a  highly  efficient  numerical  algorithm.  The  formulation  is 
sufficiently  general  so  as  to  be  applicable  to  SOS  and  bulk  device  geometries.  It 
is  also  shown  that  every  other  redistribution  problem  encountered  in  the  fabri¬ 
cation  of  VLSI  devices  is  a  special  case  of  the  field  implant  redistribution  problem. 


*This  work  was  supported  in  part  by  the  Defense  Advanced  Research  Projects  Agency  under  Contract  MDA903-80-C-0498. 
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f  Rockwell  International  Science  Center,  Thousand  Oaks,  California. 


36 


SC5271.6SA 


20  PROCESS  MOOELING  AND  SIMULATION  FOR  VLSI  DESIGN* 
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Rockwell  Intemationil  Corporation 
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SUMMARY 

As  device  geometries  are  reduced  in  site  the  use  of  two  dimensional  1201  device  and  process  simulators  as  design  aids  will  become  increasing', 
important  Two  dimensional  device  simulators  have  been  successfully  developed  during  the  past  decade  but  at  the  present  time  no  20  process 
simulator  similar  to  the  10  process  simulator,  SUPR EM  1 1  ] ,  exists.  In  a  recent  series  of  papers  12J -[41 .  the  nonlinear  boundary-value  problem 
which  governs  the  redistribution  of  implants  under  nonoxidizing  conditions,  is  solved  numerically  by  different  algorithms.  An  exact  analytical 
solution  to  this  same  problem  but  for  low  dose  implants,  has  been  obtained  by  Lee  et.  al  [51.  Also,  approximate  analytical  solutions  have  been 
obtained  [61 .  [7]  to  predict  the  redistribution  of  low  dose  implants  under  oxidizing  conditions.  In  all  these  two  dimensional  redistribution  proo 
lems  the  semiconductor  medium  is  semi-infinite  in  extent  in  the  vertical  direction  (bulk  device  type  geometry). 

The  present  trend  in  VLSI  technology  is  to  fabricate  devices  entirely  by  ion  implantations,  like  the  NMOS  enhancement  mode  device  shown  m 
Figure  1.  The  final  20  profiles  that  are  obtained  after  the  device  fabrication  schedule  it  completed,  represent  important  information  for  accurate 
electrical  device  characterization.  The  problem  of  determining  the  redistribution  of  a  field  implant  is  pertinent  since  every  other  redistribution 
problem  is  a  special  case  of  it.  Therefore,  if  an  efficient  numerical  algorithm  can  be  developed  to  solve  the  field  implant  redistribution  problem 
then  all  the  tools  which  are  necessary  to  assemble  a  SUPREM-like  2D  process  simulator,  will  be  readily  available.**  It  is  the  purpose  of  this  paper 
to  formulate,  numerically  solve,  and  apply  the  resultant  solution  to  predict  the  redistribution  of  a  field  implant.  Special  cases  pertaining  to  the 
redistribution  of  source/drain  and  enhancement  mode  channel  implants  are  aiso  considered. 

The  initial  SOS  geometry  for  defining  the  field  implant  redistribution  problem  is  shown  in  Figure  2(a).  Because  the  device  is  symmetrical  about 
the  y=0  position,  only  the  shaded  region  of  Figure  2(a)  is  used  in  formulating  the  problem.  The  window  opening  (i  yi  <  a)  in  the  impenetrable 
mask  allows  the  field  implant  to  penetrate  the  initial  thin  oxide  layer  and  distribute  itself  within  the  semiconductor  medium  in  a  two-dimensional 
manner  [81 .  When  the  photoresist  portion  of  the  impenetrable  mask  is  removed  and  the  resultant  structure  brought  into  contact  with  an  oxidizing 
environment  the  growth  of  a  field  oxide  will  take  place  with  a  lateral  penetration  of  the  bird’s  beak  under  the  mask  edge  as  shown  in  Figure  2(b) 
Applying  conservation  and  continuity  arguments  to  the  impurity  concentration  together  with  the  assumption  that  Fick's  law  191  is  applicable,  we 
derive  the  governing  boundary  value  problem  for  the  redistribution  of  a  field  implant  in  the  shaded  region  of  Figure  2(b).  For  high  dose  implants 
the  boundary  value  problem  is  highly  nonlinear  due  to  the  diffusion  coefficient's  dependence  on  impurity  concentration. 

Besides  being  nonlinear  for  high  dose  implants,  the  boundary  value  problem  is  further  complicated  by  the  fact  that  the  boundary  between  the  oxide 
and  semiconductor  media  is  moving  nonuniformly  as  a  function  of  time  and  lateral  position.  To  solve  this  boundary  value  problem  numerically  in 
a  highly  efficient  manner  the  shaded  region  in  Finite  2(b)  which  is  decreasing  in  time,  is  transformed  to  the  stationary  rectangular  shaded  region  of 
Figure  2(c).  This  mapping  was  achieved  by  a  translation-scaling  transformation.  The  stationary  rectangular  shaded  region  was  subdivided  uniformly 
in  both  the  f  and  17  directions.  On  the  resultant  computational  grid  the  classical  method  of  lines  was  used  to  discretize  the  transformed  boundary 
value  problem  into  a  stiff  system  at  coupled  ordinary  nonlinear  differential  aquations.  The  resultant  system  of  equations  was  then  solved  by  uti¬ 
lizing  the  best  currently  available  stiff  integrator,  GEARBI,  of  Hindmarsh  [101 .  This  integrator  is  a  variable  order  and  variable  step  size  method 
which  is  optimally  first  to  fifth  order  accurate  in  time.  Also,  in  GEARS),  a  successive  wer- relaxation  (SOR)  [11]  technique  is  employed,  which  is 
well  suited  to  the  present  problem  because  the  Jacobian  has  only  five  nonzero  nearest-neighbor  diagonals.  A  detail  discussion  on  the  numerical 
technique  is  presented  elsewhere  ( 121 . 

The  computer  program  based  on  the  numerical  algorithm  of  reference  (12)  was  utilized  to  predict  the  redistribution  of  an  80  keV  boron  field 
implant  of  dose  3x10*3  cm"?  during  the  growth  of  a  field  oxide  using  an  oxidizing  ambient  of  steam  at  1000°C  for  2.333  hours-  The  resultant 
equi-concentration  contours  ere  shown  in  Figure  3.  Likewise,  the  computer  program  wet  applied  to  the  special  cases  corresponding  to  the  redistri¬ 
bution  of  source/dram  and  channel  implants  end  the  results  for  the  equi-concentration  contours  are  shown  in  Figures  4  and  5,  respectively.  Figure  4 
is  for  the  redistribution  of  a  channel  profile  which  is  composed  of  a  shallow  40  keV  boron  implant  (dose  2x10**  cm'?)  superimposed  on  a  deep 
120  keV  boron  implant  (dose  1.5x10'  *  cm'?),  during  the  growth  of  a  gate  oxide  by  an  oxidizing  ambient  of  steam  at  1000°C  for  0.25  hour.  In 
Figure  5  the  redistribution  is  for  a  high  dose  5x  10*®  cm'Msource/drein)  40  k&V  arsenic  implant  which  it  being  driven-in  by  a  nonoxidizing  ambi 
ant  of  argon  et  1000°C  for  1  hour. 

The  remits  of  Figures  3,  4.  end  5  were  based  on  the  concentration  dependent  diffusion  coefficient  of  Warner  end  Wilson  [21 .  For  the  results  of 
Figure  3  the  oxide  growth  model  of  Oeel  and  Grove  (131  was  modified  empirically  to  characterize  the  lateral  position  end  time  dependence  of 
the  bird’s  beak.  Computations  were  carried  out  on  the  COC  cyber  176  computer.  To  illustrate  the  computational  efficient  of  the  numerical  algori¬ 
thm  examples  of  grid  size  end  the  resultant  computational  times  for  all  three  redistribution  problems  considered  in  this  paper  are  presented  in 
Table  I. 

'This  work  was  supported  in  part  by  the  Defense  Advanced  Research  Projects  Agency  under  Contract  MO  A903-80-C-04S8. 

•'Rockwell  International  Microelectronics  R&O  Center,  Aneheim,  California. 

f  Rockwell  International  Science  Center,  Thousand  Oaks,  California. 

f  ^Presently,  we  are  collaborating  with  the  process  modeling  group  at  Stanford  University  to  devalop  a  20  process  simulator  similar  to  SUPREM. 
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Figure  1  Qualitative  20  geometrical  sketch  of  an  NIKOS  enhancement  mode 
device  with  final  redistributed  profiles  for  rhi  field,  source/drem. 
and  channel  implants. 
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Figure  2  Initial  SOS  geometry  end  corresponding  tren stormed  geometry  met 
tekes  piece  during  me  growth  of  e  hefd  aside  ere  mown  m  (el  end 
(0).  respectively:  while  m  (cl  «  shown  me  stetionery  recfenguler 
lihaded  region)  computitionef  daman  in  tho  transformed  ((.  rj) 
coordmite  system. 


Figure  4.  Initial  end  final  ego i- boron  concentration  contain  for  the  redistribu¬ 
tion  of  e  channel  enhancement  mode  profile  ere  shown  in  la)  and  (bl. 
rupee  trvely. 
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Figure  5.  Initial  end  final  equt-ersentc  concentration  contours  lor  tho  redittr.bu- 
tion  of  t  source/drein  implont  are  shown  in  (a)  end  (bl.  rtspecnvfiy 
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